
version 14.2
clear all
set more off

*** addos needed:
*cap ssc install spmap 
*cap ssc install shp2dta
*cap ssc install mif2dta 


* Gemeinde
shp2dta using "data/orig/vg250_gem.shp", database("data/ao_gem_map.dta") coordinates("data/orig/GER_ao_gem_coord.dta") genid(id) gencentroids(centroid) replace 
use "data/ao_gem_map.dta" , clear
duplicates tag AGS, gen(duplo)
sort AGS EWZ
bys AGS (EWZ): drop if _n==1 & duplo==1
drop duplo
ren AGS ao_gem
destring ao_gem , replace
format ao_gem %12.0g
sort ao_gem
save "data/ao_gem_map.dta" , replace


/* GEMEINDE: input coordinates of border crossings and calculate distances*/
use "data/ao_gem_map.dta", replace

	/* Transform coordinates into WGS84 */ 
	* NOTE: Code adopted from DELPHI code from http://www.delphi-treff.de/tipps/mathematik/einheiten/geographische-in-gauss-krueger-koordinaten-umrechnen/
 
  local rho= 180 / _pi

  local e2 = 0.0067192188

  local c = 6398786.849
 
  gen double mKen = int(x_centroid / 1000000)

  gen double rm = x_centroid - mKen * 1000000 - 500000
 
  gen double bI = y_centroid / 10000855.7646

  gen double bII = bI * bI

  gen double bf = 325632.08677 *bI *((((((0.00000562025   * bII - 0.00004363980)  * bII + 0.00022976983)  * bII - 0.00113566119)  * bII + 0.00424914906)  * bII - 0.00831729565)  * bII + 1)

  replace bf = bf / 3600 / `rho'

  gen double co = cos(bf)

  gen double g2 = `e2' * (co * co)

  gen double g1 = `c' / sqrt(1 + g2)

  gen double t = sin(bf) / cos(bf)

  gen double fa = rm / g1

  gen double gb = bf        - fa * fa * t * (1 + g2) / 2  + fa * fa * fa * fa * t * (5 + 3 * t * t + 6 * g2 - 6 * g2 * t * t) / 24

  replace gb = gb * `rho'

  gen double dl = fa        - fa * fa * fa * (1 + 2 * t * t + g2) / 6 + fa * fa * fa * fa * fa * (1 + 28 * t * t + 24 * t * t * t * t) / 120

  gen double gl = dl *`rho' / co + mKen * 3

  ren gl gl_centroid
  ren gb gb_centroid

	label var x_centroid "GK3 Geocoord. East"
	label var y_centroid "GK3 Geocoord. North"
	label var gb_centroid "WGS84 Latitude"
	label var gl_centroid "WGS84 Longitude"
  
  	* Clean
  	drop mKen - fa 
  	drop dl

	* input coordinates of German-Czech border crossings

input x_CROSSING y_CROSSING gb_CROSSING gl_CROSSING
	3725543.70 5580975.44 50.321573 12.166271
	3726377.23 5568884.62 50.212724 12.170734
	3733041.85 5555103.61 50.086446 12.255535
	3739195.17 5550145.18 50.039502 12.338276
	3755009.80 5535199.56 49.898836 12.548733
	3749366.76 5525241.55 49.811847 12.463932
	3754608.83 5508003.70 49.654960 12.525387
	3769566.16 5482304.78 49.417897 12.714558	
	3779865.53 5474226.57 49.340783 12.850513
	3788735.30 5473913.58 49.333848 12.972050
	3806954.62 5451351.30 49.122646 13.204823
	3845097.75 5426586.16 48.880522 13.704700
end
	label var x_CROSSING "GK3 Geocoord. East of Border Crossings"
	label var y_CROSSING "GK3 Geocoord. North of Border Crossings"
	label var gb_CROSSING "WGS84 Latitude of Border Crossings"
	label var gl_CROSSING "WGS84 Longitude of Border Crossings"
	
	* input coordinats of Plzen (major city in West Czech Republic)
input x_Plzen y_Plzen gb_Plzen gl_Plzen
	3815539.56 5521479.06 49.747029 13.377821 
end	

	label var x_Plzen "GK3 Geocoord. East of Plzen/Pilsen"
	label var y_Plzen "GK3 Geocoord. North of Plzen/Pilsen"
	label var gb_Plzen "WGS84 Latitude of Plzen/Pilsen"
	label var gl_Plzen "WGS84 Longitude of Plzen/Pilsen"
	
	* install nearstat if not yet installed and calculate minimum distance to tagged border-municipalities, to Plzen
	mata: mata mlib index

	nearstat gb_centroid gl_centroid, near(gb_CROSSING gl_CROSSING) distvar(distance) kth(1)
	nearstat gb_centroid gl_centroid, near(gb_Plzen gl_Plzen) distvar(distance_Plzen) kth(1)

	label var distance "distance (airline) to nearest border crossing (in km), accuracy munic. level"
	label var distance_Plzen "distance (air line) to Plzen, accuracy munic. level"
sort ao_gem
rename ao_gem  ao_gem_imp

save "data/ao_gem_map.dta", replace

exit
